Supplementary Fig 8 - GEP usage per thymic lineage
Salomé Carcy
29 July 2024
1 IMPORT
1.1 Import librairies
library(ggplot2)
library(cowplot)
library(tidyverse)
library(dplyr)
library(Seurat)
setwd("~/Projects/HumanThymusProject/")
source("~/Projects/HumanThymusProject/scripts/colors_universal.R")1.2 Import data
# import seurat objects
seur.nkt <- readRDS("./data_github/seurat_objects_thymus/seurat_thymus_inkt.rds")
seur.mait <- readRDS("./data_github/seurat_objects_thymus/seurat_thymus_mait.rds")
seur.gdt <- readRDS("./data_github/seurat_objects_thymus/seurat_thymus_gdt.rds")
seur.cd4 <- readRDS("./data_github/seurat_objects_thymus/seurat_thymus_cd4.rds")
seur.cd8 <- readRDS("./data_github/seurat_objects_thymus/seurat_thymus_cd8.rds")
# sanity check seurat objects correspond to Fig 2
DimPlot(seur.nkt, group.by="clusters_per_lineage", label=T, cols=cols_thym_nkt)
DimPlot(seur.mait, group.by="clusters_per_lineage", label=T, cols=cols_thym_mait)
DimPlot(seur.gdt, group.by="clusters_per_lineage", label=T, cols=cols_thym_gdt)
# import GEP usage dataframe
gep_usage <- read.table("./data_github/cNMF/non_imputed_cNMF_allcells.usages.k_12.dt_0_02.consensus.txt", header=T)
colnames(gep_usage) <- paste0("GEP", c(2,5,3,1,4,12,6,7,8,10,9,11), "_usage")
gep_usage <- gep_usage[,paste0("GEP", 1:12, "_usage")]
head(gep_usage)## GEP1_usage GEP2_usage GEP3_usage GEP4_usage GEP5_usage
## CD4_1_Thymus_878524 0.0000000 0.00000000 0.01697069 0.02701780 0.003293012
## CD4_1_Thymus_679405 0.0000000 0.00000000 0.00000000 0.00000000 0.000000000
## CD4_1_Thymus_281188 0.0000000 0.01499587 0.02538062 0.01897941 0.000000000
## CD4_1_Thymus_285601 0.3893558 0.19701648 0.01374016 0.03301096 0.009603295
## CD4_1_Thymus_174537 0.0000000 0.00000000 0.01813514 0.02750585 0.001959355
## CD4_1_Thymus_128406 0.0000000 0.02029776 0.02094953 0.01324090 0.001427132
## GEP6_usage GEP7_usage GEP8_usage GEP9_usage GEP10_usage
## CD4_1_Thymus_878524 0 0.0000000000 0.01670451 0.06574529 0.14537527
## CD4_1_Thymus_679405 0 0.0000000000 0.00000000 0.00000000 0.00000000
## CD4_1_Thymus_281188 0 0.0002823688 0.01075830 0.09749147 0.21138272
## CD4_1_Thymus_285601 0 0.0178349300 0.22377394 0.05961335 0.01831705
## CD4_1_Thymus_174537 0 0.0000000000 0.01610514 0.06646363 0.15172973
## CD4_1_Thymus_128406 0 0.0055636405 0.01326269 0.08876818 0.21684438
## GEP11_usage GEP12_usage
## CD4_1_Thymus_878524 0.655783300 0.00000000
## CD4_1_Thymus_679405 1.147190900 0.00000000
## CD4_1_Thymus_281188 0.563456800 0.00000000
## CD4_1_Thymus_285601 0.008145613 0.03549316
## CD4_1_Thymus_174537 0.649509700 0.00000000
## CD4_1_Thymus_128406 0.565640400 0.000000002 ANALYSIS
2.1 Add GEP1-GEP11 usage to seurat objects
# check that rownames are in same order
# table(rownames(gep_usage[rownames(gep_usage) %in% rownames(seur.cd4@meta.data),])==rownames(seur.cd4@meta.data), useNA="ifany")
# table(rownames(gep_usage[rownames(gep_usage) %in% rownames(seur.cd8@meta.data),])==rownames(seur.cd8@meta.data), useNA="ifany")
# table(rownames(gep_usage[rownames(gep_usage) %in% rownames(seur.nkt@meta.data),])==rownames(seur.nkt@meta.data), useNA="ifany")
# table(rownames(gep_usage[rownames(gep_usage) %in% rownames(seur.mait@meta.data),])==rownames(seur.mait@meta.data), useNA="ifany")
# table(rownames(gep_usage[rownames(gep_usage) %in% rownames(seur.gdt@meta.data),])==rownames(seur.gdt@meta.data), useNA="ifany") # all TRUE
# add GEP usage to seurat objects
gep_colnames <- paste0("GEP", 1:11, "_usage")
seur.cd4@meta.data[,gep_colnames] <- gep_usage[rownames(gep_usage) %in% rownames(seur.cd4@meta.data), gep_colnames]
seur.cd8@meta.data[,gep_colnames] <- gep_usage[rownames(gep_usage) %in% rownames(seur.cd8@meta.data), gep_colnames]
seur.nkt@meta.data[,gep_colnames] <- gep_usage[rownames(gep_usage) %in% rownames(seur.nkt@meta.data), gep_colnames]
seur.mait@meta.data[,gep_colnames] <- gep_usage[rownames(gep_usage) %in% rownames(seur.mait@meta.data), gep_colnames]
seur.gdt@meta.data[,gep_colnames] <- gep_usage[rownames(gep_usage) %in% rownames(seur.gdt@meta.data), gep_colnames]2.2 Plot GEP usage per thymic lineage
# Plot
seur_vector <- list("iNKT"=seur.nkt, "MAIT"=seur.mait, "GDT"=seur.gdt, "CD4"=seur.cd4, "CD8"=seur.cd8)
plist_massive <- list()
for(gep in gep_colnames){
print(gep)
# get row of plots (for one GEP)
plist <- list()
for(i_seur in names(seur_vector)){
print(i_seur)
seur <- seur_vector[[i_seur]]
# get legend once
if(gep=="GEP1_usage" & i_seur=="iNKT"){
plegend <- ggpubr::get_legend(
SCpubr::do_FeaturePlot(seur, features=gep, ncol=1, legend.position="right", legend.title="GEP usage", border.size=1, pt.size=2, order=T)+
scale_color_viridis_c(limits=c(0,1.25), option="B")+
theme(plot.background = element_rect(color = "black"))
)
}
# make plot
legendpos <- "none"
plt_title <- ""
if(gep=="GEP1_usage"){plt_title=i_seur}
p <- SCpubr::do_FeaturePlot(seur, features=gep, ncol=1, legend.position=legendpos, border.size=1, pt.size=2, order=T, plot.title=plt_title)+
scale_color_viridis_c(limits=c(0,1.25), option="B")+
theme(plot.background = element_rect(color = "black", linewidth = 2),
plot.margin=unit(c(0,5.5,0,5.5), "points"))
plist[[i_seur]] <- ggrastr::rasterise(p, layers="Point", dpi=300)
}
prow_gep <- plot_grid(plotlist=plist, nrow=1, scale=0.9)
plist_massive[[gep]] <- prow_gep
}## [1] "GEP1_usage"
## [1] "iNKT"
## [1] "MAIT"
## [1] "GDT"
## [1] "CD4"
## [1] "CD8"
## [1] "GEP2_usage"
## [1] "iNKT"
## [1] "MAIT"
## [1] "GDT"
## [1] "CD4"
## [1] "CD8"
## [1] "GEP3_usage"
## [1] "iNKT"
## [1] "MAIT"
## [1] "GDT"
## [1] "CD4"
## [1] "CD8"
## [1] "GEP4_usage"
## [1] "iNKT"
## [1] "MAIT"
## [1] "GDT"
## [1] "CD4"
## [1] "CD8"
## [1] "GEP5_usage"
## [1] "iNKT"
## [1] "MAIT"
## [1] "GDT"
## [1] "CD4"
## [1] "CD8"
## [1] "GEP6_usage"
## [1] "iNKT"
## [1] "MAIT"
## [1] "GDT"
## [1] "CD4"
## [1] "CD8"
## [1] "GEP7_usage"
## [1] "iNKT"
## [1] "MAIT"
## [1] "GDT"
## [1] "CD4"
## [1] "CD8"
## [1] "GEP8_usage"
## [1] "iNKT"
## [1] "MAIT"
## [1] "GDT"
## [1] "CD4"
## [1] "CD8"
## [1] "GEP9_usage"
## [1] "iNKT"
## [1] "MAIT"
## [1] "GDT"
## [1] "CD4"
## [1] "CD8"
## [1] "GEP10_usage"
## [1] "iNKT"
## [1] "MAIT"
## [1] "GDT"
## [1] "CD4"
## [1] "CD8"
## [1] "GEP11_usage"
## [1] "iNKT"
## [1] "MAIT"
## [1] "GDT"
## [1] "CD4"
## [1] "CD8"plot_grid(
plot_grid(plotlist=plist_massive, ncol=1, align="h"),
ggpubr::as_ggplot(plegend), ncol=2, rel_widths = c(5, .5), scale=c(0.95, 5)
)
3 SESSION INFO
sessionInfo()## R version 4.3.1 (2023-06-16)
## Platform: x86_64-apple-darwin20 (64-bit)
## Running under: macOS Sonoma 14.2.1
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/4.3-x86_64/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.3-x86_64/Resources/lib/libRlapack.dylib; LAPACK version 3.11.0
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## time zone: Europe/Paris
## tzcode source: internal
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] Seurat_5.0.1 SeuratObject_5.0.1 sp_2.1-2 lubridate_1.9.3
## [5] forcats_1.0.0 stringr_1.5.1 dplyr_1.1.4 purrr_1.0.2
## [9] readr_2.1.5 tidyr_1.3.0 tibble_3.2.1 tidyverse_2.0.0
## [13] cowplot_1.1.2 ggplot2_3.4.4
##
## loaded via a namespace (and not attached):
## [1] RColorBrewer_1.1-3 rstudioapi_0.15.0 jsonlite_1.8.8
## [4] magrittr_2.0.3 ggbeeswarm_0.7.2 spatstat.utils_3.0-4
## [7] farver_2.1.1 rmarkdown_2.25 vctrs_0.6.5
## [10] ROCR_1.0-11 Cairo_1.6-2 spatstat.explore_3.2-5
## [13] rstatix_0.7.2 htmltools_0.5.7 broom_1.0.5
## [16] sass_0.4.8 sctransform_0.4.1 parallelly_1.36.0
## [19] KernSmooth_2.23-22 bslib_0.6.1 htmlwidgets_1.6.4
## [22] ica_1.0-3 plyr_1.8.9 plotly_4.10.3
## [25] zoo_1.8-12 cachem_1.0.8 igraph_1.6.0
## [28] mime_0.12 lifecycle_1.0.4 pkgconfig_2.0.3
## [31] Matrix_1.6-5 R6_2.5.1 fastmap_1.1.1
## [34] fitdistrplus_1.1-11 future_1.33.1 shiny_1.8.0
## [37] digest_0.6.34 colorspace_2.1-0 patchwork_1.2.0
## [40] tensor_1.5 RSpectra_0.16-1 irlba_2.3.5.1
## [43] ggpubr_0.6.0 labeling_0.4.3 progressr_0.14.0
## [46] fansi_1.0.6 spatstat.sparse_3.0-3 timechange_0.2.0
## [49] httr_1.4.7 polyclip_1.10-6 abind_1.4-5
## [52] compiler_4.3.1 withr_2.5.2 backports_1.4.1
## [55] viridis_0.6.4 carData_3.0-5 fastDummies_1.7.3
## [58] highr_0.10 ggsignif_0.6.4 MASS_7.3-60
## [61] tools_4.3.1 vipor_0.4.7 lmtest_0.9-40
## [64] beeswarm_0.4.0 httpuv_1.6.13 future.apply_1.11.1
## [67] goftest_1.2-3 glue_1.7.0 nlme_3.1-164
## [70] promises_1.2.1 grid_4.3.1 Rtsne_0.17
## [73] cluster_2.1.6 reshape2_1.4.4 generics_0.1.3
## [76] gtable_0.3.4 spatstat.data_3.0-3 tzdb_0.4.0
## [79] data.table_1.14.10 hms_1.1.3 car_3.1-2
## [82] utf8_1.2.4 spatstat.geom_3.2-7 RcppAnnoy_0.0.21
## [85] ggrepel_0.9.5 RANN_2.6.1 pillar_1.9.0
## [88] spam_2.10-0 RcppHNSW_0.5.0 later_1.3.2
## [91] splines_4.3.1 lattice_0.22-5 survival_3.5-7
## [94] deldir_2.0-2 tidyselect_1.2.0 miniUI_0.1.1.1
## [97] pbapply_1.7-2 knitr_1.45 gridExtra_2.3
## [100] SCpubr_2.0.2 scattermore_1.2 xfun_0.41
## [103] matrixStats_1.2.0 stringi_1.8.3 lazyeval_0.2.2
## [106] yaml_2.3.8 evaluate_0.23 codetools_0.2-19
## [109] cli_3.6.2 uwot_0.1.16 xtable_1.8-4
## [112] reticulate_1.34.0 munsell_0.5.0 jquerylib_0.1.4
## [115] Rcpp_1.0.12 globals_0.16.2 spatstat.random_3.2-2
## [118] png_0.1-8 ggrastr_1.0.2 parallel_4.3.1
## [121] ellipsis_0.3.2 assertthat_0.2.1 dotCall64_1.1-1
## [124] listenv_0.9.0 viridisLite_0.4.2 scales_1.3.0
## [127] ggridges_0.5.5 leiden_0.4.3.1 rlang_1.1.3